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Abstract 

Tlie analysis of gene networl< robustness to noise and mutation is important for fundamental and practical reasons. 
Robustness refers to the stability of the equilibrium expression state of a gene network to variations of the initial 
expression state and network topology. Numerical simulation of these variations is commonly used for the 
assessment of robustness. Since there exists a great number of possible gene network topologies and initial states, 
even millions of simulations may be still too small to give reliable results. When the initial and equilibrium expression 
states are restricted to being saturated (i.e., their elements can only take values 1 or — 1 corresponding to maximum 
activation and maximum repression of genes), an analytical gene network robustness assessment is possible. We 
present this analytical treatment based on determination of the saturated fixed point attractors for sigmoidal function 
models. The analysis can determine (a) for a given network, which and how many saturated equilibrium states exist 
and which and how many saturated initial states converge to each of these saturated equilibrium states and (b) for a 
given saturated equilibrium state or a given pair of saturated equilibrium and initial states, which and how many gene 
networks, referred to as viable, share this saturated equilibrium state or the pair of saturated equilibrium and initial 
states. We also show that the viable networks sharing a given saturated equilibrium state must follow certain patterns. 
These capabilities of the analytical treatment make it possible to properly define and accurately determine robustness 
to noise and mutation for gene networks. Previous network research conclusions drawn from performing millions of 
simulations follow directly from the results of our analytical treatment. Furthermore, the analytical results provide 
criteria for the identification of model validity and suggest modified models of gene network dynamics. The yeast 
cell-cycle network is used as an illustration of the practical application of this analytical treatment. 

Keywords: Robustness of gene networks; Fixed point attractor; Sigmoidal function; Yeast cell-cycle network 



1 Introduction 

A subset of genes in a cell whose protein products mutu- 
ally regulate one another's expression at the transcrip- 
tional level will be referred to as a 'gene network' in this 
paper. The concentration of proteins encoded by the genes 
changes in time due to auto- and cross-regulation by the 
gene products. Each of such network is considered as a 
dynamical system. Gene networks must be robust with 
respect to ever-changing environments. Robustness here 
refers to the ability of a gene network to respond to short- 
term changes in the environment and quickly return to 
its functional steady state. Moreover, a gene network itself 
may endure small structural changes and mutations, while 
still retaining its desired steady state. The robustness of 
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gene networks depends on their topology with some net- 
works being more stable than others. The analysis of 
the relationship between the topology of a gene network 
and its robustness to environmental and structural per- 
turbations is important both theoretically and practically 
[1-10]. 

Recently, Wagner and coworkers considered the robust- 
ness of gene networks [11-14], and a similar assessment 
was given by Cho et al. [15] In these works, a simpli- 
fied model proposed by Wagner was used to describe the 
dynamics of the gene expression states [11,12]. Let G = 
(Gi, . . . , G„) represent the n genes in a network. The con- 
centration of proteins encoded by the genes (Gi, . . . , G„) 
is denoted by P = (Pi, . . . ,Pn)- For computational con- 
venience, the admissible concentration range for each Pi 
is normalized and restricted to the interval [0, 1], where 
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Pi = 1 corresponds to the maximum possible concen- 
tration, i.e., the corresponding gene Gi is in a state of 
maximum transcriptional activation. It is also assumed 
that Pi = 0.5 means that gene i is 50% 'on'. The dynamics 
of the expression states of the genes in a network is often 
described by some sigmoidal function j^cW 

Piit + r}=gc^^WijPjit}j, 0 = 1,2,...,«) (1) 

where r is a time constant characteristic of the process 
under consideration. In some work, r was set to be 1. 
The constant Wij e 9t describes the strength of interac- 
tion (i.e., transcriptional regulation) of the product of gene 
j with gene i, i.e., the degree of transcriptional activation 
{wij > 0), repression (tVy < 0), or absence [wij = 0). These 
constants define a matrix of connectivities W = (wi/) 
within the network. To facilitate the analytical treatment, 
the variable transformation 



S = 2P- (1 



1)^ 



(2) 



is employed. Using the sigmoidal function a proposed by 
Siegal [16] and Cho [15] for^c. (1) becomes 

Siit + r) = a ^^WijSj(t)j 

^ 1, (i = 1,2,...,m) 





n 


1 + exp 









(3) 

with Si g[ — 1, 1], and 5, = 0 corresponding to 50% of gene 
/ being 'on'. Notwithstanding the simplicity of (3), vari- 
ants of this model have been successfully used to study (a) 
the robustness of gene regulatory networks [12,16,17], (b) 
the role of robustness in evolutionary innovation [18,19], 
and (c) how recombination can produce negative epistasis 
[20]. 

Mjolsness at al. [21] proposed a model 



' dt 



--Sa\J2T''^yl+h'']-XaVf, (« = 1,2,...,«) 



(4) 



to describe the dynamics for each element of primitive 
objects V/ (cells, nuclei, fibers, and synapses), where ga is 
a sigmoidal threshold function, T"^ is similar to in (3); 
v'f, denote the elements of vector v,; h" determines the 
threshold of The long-time behavior of this system has 
been studied, and, in some cases, is controlled by a simple 
limit set 



(5) 



which is similar to (3) with the additional parameters 
and W^. This model has been successfully applied to treat 
the blastoderm of Drosophila melanogaster [21]. 

Similarly, Mendoza and Alvarez-Buylla [22] used a 
model 



where H is the Heaviside step function 
H{x) = 



1, if;« > 0, 
0, iix < 0 



(6) 



(7) 



to describe the dynamics of a genetic regulatory net- 
work for Arabidopsis thaliana flower morphogenesis. 
This model is also similar to (3) except that the sigmoidal 
function a is replaced by the Heaviside step function and 
a threshold parameter 0i is included. 

All these models present simplified descriptions of 
gene network dynamics. Nevertheless, the models are 
still useful for obtaining insights into the dynamics of 
gene networks. In the following analysis for gene network 
robustness, we will employ the sigmoidal function model 
in (3), and its modification with threshold parameters. 

The robustness of a gene network specified by W to 
noise (environmental) and mutation (structural) pertur- 
bations may be expressed as the stability of the final 
equilibrium (or steady) expression state S(oo) obtained 
from the solution of (3), respectively, to changes of ini- 
tial expression state S(0) and to changes of W. A complete 
and reliable robustness analysis would seem to call for 
an exhaustive sampling over the space of possible ini- 
tial states for a given network and then repeating the 
same simulations for all possible networks. This task is 
infeasible as there are many possible networks W, and 
each W may have many initial/final equilibrium expres- 
sion states. Consider a simple case where S, (0) and S/(oo) 
{i = 1, 2, . . .,«) can only take values —1,1 and Wij can 
only take values —1, 0, 1. In this case, there are 2" possible 

2 

initial/final states and 3" possible gene networks. For a 
modest network of size n = 20, there are 2^° = 1, 048, 576 
initial/final expression states and 3'**^^ ^ 7 x 10^^^ gene 
networks. Even if one arbitrarily makes the restriction that 
75% of the Wij interactions are zero, and that the remain- 
ing 25% of the Wij interactions can only take nonzero 
values —1,1 to reduce the possible number of networks 
[13], for « = 20 there are still C^^ x 2^°° possible gene net- 
works. Further restriction may also be applied to reduce 
the possible number of initial expression states [13]. Even 
under these restrictions, solving (3) for all possible initial 
expression states and gene networks is still infeasible, and 
only a small fraction of them can be randomly sampled 
and simulated. Such limitations leave open the reliability 
of the conclusions obtained from the simulations. 
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Previous work [13] concerned networks with connec- 
tivity W whose expression dynamics start from a pre- 
specified initial state S(0) at some time t = 0 and arrive 
at a prespecified stable equilibrium or 'target' expression 
state S(oo); these networks are referred to as 'viable'. 
Then, the values of some elements of S(0) or W were 
changed for each viable network to check whether S(oo) 
is reached. These studies entailed performing millions of 
simulations with different network topologies and initial 
expression states. Although the number of simulations is 

2 

much smaller than 2" and 3" for n = 20, some specific 
conclusions were obtained. First, the fraction of viable 
networks, that is, networks that arrive at a prespecified 
target expression state S(oo) given an initial gene expres- 
sion state S(0) to the total number of possible networks, 
is generally very small. For moderately sized networks of 
n = 20 genes (with the number M of nonzero w,y set to 
be 200, and the fraction d of elements different between 
S(0) and S(oo) set to 0.5), the fraction of viable networks 
was found to be Vf = 5.1 x 10"^ ± 1.7 x 10"^°. Due 

2 

to the large numbers 2" and 3" , the qualitative correct- 
ness of this conclusion is clear. Since there are 2" possible 
equilibrium states, even if we only consider the factor of 
expression states, the probability that a network W arrives 
at a prespecified S(oo) is expected on the order of 1/2". 
The viable networks in this prior work could be organized 
as a graph with each node corresponding to a network 
of a given topology, and two nodes are connected by an 
edge if they differ by a single regulatory interaction (i.e., 
they differ in one element Wy). Remarkably, this graph is 
connected and can be easily traversed by gradual changes 
of the network topology. Thus, highly robust topologies 
can evolve from topologies with low robustness through 
gradual Darwinian topological changes. These results are 
claimed to be valid for discrete and continuous taking 
values [—1,0,1] and over the interval [—a, a], respec- 
tively. While simulations are valuable, they do not provide 
a complete picture. 

Ciliberti et al. [13] considered the case where each 
element of the initial and equilibrium expression states, 
S(0) and S(oo), can only take the values 1 or —1 cor- 
responding to maximum and minimum possible protein 
concentrations. We call them saturated expression states. 
Under this condition, the present paper provides an ana- 
lytical robustness assessment of gene networks whose 
dynamics can be described by (3) and its modification 
with threshold parameters. This analysis can determine 
(a) for a given network, which and how many saturated 
equilibrium states exist, and which and how many sat- 
urated initial states converge to each of these saturated 
equilibrium states; (b) for a given saturated equilibrium 
state, or a given pair of saturated equilibrium and initial 
states, which and how many gene networks, referred to as 
viable, share this saturated equilibrium state or the pair 



of saturated equilibrium and initial states. We also show 
that the viable networks sharing a given saturated equi- 
librium state must follow certain patterns. These capa- 
bilities of the analytical treatment make it possible to 
properly define and accurately determine robustness to 
noise and mutation for gene networks. Previous network 
research conclusions drawn from performing millions of 
simulations follow directly from the results of our analyti- 
cal treatment. Furthermore, the analytical results provide 
criteria for identification of model validity and suggest 
modified models of gene network dynamics. 

The paper is organized as follows: Section 2 first defines 
the saturated state and saturated fixed point attractor for 
dynamics (3), and then gives the necessary and sufficient 
condition for a gene network to have a given saturated 
equilibrium state. Sections 3 and 4 analyze the robustness 
to noise and mutation, respectively. Section 5 proposes a 
modification of dynamics (3) with threshold parameters. 
Section 6 gives an illustration of the practical applica- 
tion of this analytical treatment: the model construction 
of the yeast cell-cycle network and its robustness assess- 
ment. The details of the treatment of the yeast cell-cycle 
network are given in an Additional file 1: Supplemen- 
tary information. Finally, Section 7 presents conclusions. 
Mathematical proofs of the theorems in the main text are 
given in the Appendix. 

2 Saturated states and fixed point attractors 

In this work, the initial and equilibrium expression states 

S,(0) and Si{oo) for gene i can only be either active {Si{0) 
and Siioo) = 1) or inactive (5,(0) and SKoo) = — 1) 
[11,12,14]. The initial and equilibrium expression states 
with S,(0) and Si(oo) = ±1 are referred to as saturated 
initial and equilibrium expression states. Under the condi- 
tion that the initial and equilibrium expression states 5(0) 
and S(oo) are saturated, we may analyze the robustness 
and evolvability of gene networks analytically, as explained 
below. 

2.1 Saturated sigmoidal function 

A continuous function/(x) defined on satisfying 

(fl) f{x) = \^x> l,f(x) = -1 iix < -1, 

(f2) f{x) is a strictly increasing and continuous function 

for;*; €[-1,1] and/(0) = 0, 

is called a saturated sigmoidal function. Furthermore, if 

(f3) fix) >xforx& (0, 1] and fix)<xforx&[-l, 0), 

we call/ (x)] a dissipative saturated sigmoidal function. 

Note that for the particular sigmoidal function cr^ (x) in 
domain [-1, 1] 

x(t + T) = a^ixit)) = ^^^i^ - 1, (8) 
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the conditions (fl) to (f3) are approximately satisfied 
when jS is sufficiently large. For example, when /J = 5 and 
10, we have apil) = 0.9866 and 0.9999 (approximately 1); 
and ap(-l) = -0.9866 and -0.9999 (approximately -1), 
respectively. Therefore, in numerical simulation, ap(x) 
can be considered as a dissipative saturated sigmoidal 
function for a sufficiently large p. We refer to /J = 5 and 10 
as having 0.99 and 0.9999 confidence levels for the dissi- 
pative saturated sigmoidal function (8), because |cr|6(±l)| 
is equal to 0.99 and 0.9999, respectively (see Figure 1). In 
the sequel, we set /J > 5. 



It is easy to see that a saturated state S is a saturated 
fixed point of dynamics (3) with 0.99 (jS = 5) or 0.9999 
(/8 = 10) confidence level when 



^ Wij sign(5,) > P, if Si = 1, 



/=i 



^ Wij signiSj) < -p, if Si = -1, 
7=1 



(9) 



(10) 



and 



2.2 Necessary and sufficient condition for a saturated 
state S to be an equilibrium state or a fixed point 
attractor of dynamics (3) with a given W 

When the equilibrium states S(oo) are saturated, the anal- 
ysis of robustness and evolvabUity for gene networks can 
be readily performed by utilizing Feng and Tirozzi's treat- 
ment for neural networks [23]. 

Definition 1. A saturated state S in [—1,1]" is called 
a saturated fixed point attractor (or saturated equilib- 
rium state) of dynamics (3), if there exists a nonempty 
neighborhood B(S) of S such that 



lim S(t)=S 

t— >CX) 



for S(0) e B(S) and ^jLi WijSj ^ 0, for all i. 



lim Si(t) 



1, if Si = 1 
-1, if Si = -1 



(i= 1,2,..., n), 



(11) 



i.e., S is a saturated fixed point with 0.99 [p = 5) or 0.9999 
[p = 10) confidence level. When 



^ Wij sign (5,) 



< 5, 



for any /, then —0.99 < Si{t + r) < 0.99 for any t, and the 
gene network cannot have a saturated fixed point. 
Let S be a saturated state. We define 



r(S) = {i. Si = 1}, /-(S) = {i. Si = -1} 



(12) 



which denote, respectively, the two sets of all integers in 
{1, 2, ... , n] with Si = 1 and 5, = -1. 
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Theorem 1. The necessary and sufficient condition for 
a saturated state S to be an equilibrium expression state or 
a fixed point attractor of dynamics (3) with a given matrix 



Wis 






E 


;€/+(S) 


;e/-(S) 




E 




/e/-(S) 


or equivalently 





J2 Wii>P, if«e/+(S), 



(13) 



5; I ^M/yS,- j >/3, 0 = 1,2,...,«). 



(15) 



Proof. If a saturated state S is a fixed point of (3), it 
must satisfy (9,10), i.e., (13,14). If (13,14) are satisfied, so 
are (9,10), then (11) is satisfied and S is a fixed point. A 
saturated state S satisfying (13,14) (or equivalently (15)) is 
a fixed point attractor. The proof is given in Theorem A2 
in the Appendix . □ 

There are a total 2" saturated states S. Since (15) only 
involves simple multiplication and summation, for a mod- 
est n one can test all 2" saturated vectors S for a given W 
and find all of its saturated fixed point attractors without 
iteratively solving the sigmoidal function (3). For « = 11 
and 2" = 2, 048, the test takes 0.01 s by Matlab on a Dell 
Precision Workstation T3400. 

Example 1. Consider the network model proposed by 
Azevedo et al. [20] for the gap gene system of Drosophila 
melanogaster shown in Figure 2. 

The authors obtained the equilibrium state Si = 
(—11 — 11) by iteratively solving a sigmoidal function sim- 
ilar to (3). Using (9,10), we can determine this solution by 
noting 



W-Si 



2 


-1 


-1 


-5' 




"-1' 




' -7 


8 


1 


0 


I 




1 




10 


1 


-1 


1 


-4 




-1 




-5 


0 


-1 


-6 


1 




1 




6 



which shows that the necessary and sufficient condition is 
satisfied with 0.99 confidence level 



J^'^iA^^' « = 2,4e/+(Si), 



7=1 

4 

WijSj < -5, J = 1, 3 e /- (Si). 

There are 2'* = 16 saturated states. Similar tests for 
the other 15 saturated states were performed. The case 
S2 = —Si = (1 — 11 — 1) is the only other saturated 
equilibrium state for the network. 



3 Robustness to noise 

Robustness to noise may be assessed for (a) each satu- 
rated equilibrium expression state or (b) a specified pair 
of saturated equilibrium and initial expression states. In 
case 1, we need to compare how many of 2" possible 
saturated initial expression states converge to each sat- 
urated equilibrium expression state; in case 2, we need 
to determine how many neighbours (differing only in one 
element) of the S(0) converge to the same saturated equi- 
librium expression state for a given W [13]. In either case, 
we need to establish the condition under which a satu- 
rated initial expression state converges to a given saturated 
equilibrium expression state of W. 

3.1 Relationship between saturated Initial expression 
states and saturated equilibrium expression states 
Theorem 2. If S is a saturated equilibrium expression 
state (or a fixed point attractor) of dynamics (3) with a 
given W, so is — S. 

Proof. Since S is a saturated fixed point attractor of (3) 
with a given W, then (15) is satisfied. Multiplying by — 1 
within and outside the parentheses in (15) will not change 
its right-hand side: 

(-1)5; l^(-l) ^ M/yS; j >p, a =1,2,.... n), 
i-Si) (-5,0 j >p, a =1,2,..., n) 

(16) 

which proves that — S is also a saturated fixed point attrac- 
tor, i.e., a saturated equilibrium expression state for W. 
Example 1 demonstrates its validity. □ 

Corollary 1. Dynamics (3) either does not have a satu- 
rated equilibrium expression state, or has an even number 
of saturated equilibrium expression states. 

This result can be obtained immediately from Theorem 2. 

Example 2. Consider a gene network given by 



W ■■ 



0 


1 


1 


-1 


-1 


-1 


1 


0 


1 


-1 


-1 


-1 


1 


1 


0 


-1 


-1 


-1 


-1 


-1 


-1 


0 


1 


1 


-1 


-1 


-1 


1 


0 


1 


-1 


-1 


-1 


1 


1 


0 



(17) 



with the discrete values [-1, 0, 1] for Wij, and there is no 
auto-regulation (wu = 0). The six genes are separated 
into two groups {1, 2, 3} and {4, 5, 6}. The regulations are 
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kni 



Kr 



Time t 



hb 



kni 



Kr 




kni 



Kr 




t+^ 



t + 2 



Figure 2 The network model for the gap gene system of Drosophila melanogaster (reproduced from [20]). (a) Network representation of the 
regulatory interactions between four gap genes (gt, giant; hb, hunchbacl<; kni, l<nirps; Kr, Kruppel). (b) Corresponding matrix W. (c) graphic 
representation of the gene expression states of each gap gene over three successive time steps of a sigmoidal function similar to (3) where a filled 
box denotes 'on' (1), and an open box denotes 'off (—1). 



activating within each group, but repressing between the 
two groups. For such a simple system, it is easy to find 
by observation that amongst the 2^ = 64 saturated states 
only two states 

51 = (1 1 1 -1 -1 -1), 

52 = (-1 -1 -1 1 1 1) 

satisfy the necessary and sufficient condition (15) to be 
saturated equilibrium expression states with 0.99 confi- 
dence level. An examination of (15) for all 64 saturated 
states proved this to be the case. Moreover, 

S2 = -Si 
satisfies Theorem 2 and Corollary 1. 

Theorem 3. If S(t) converges to a saturated equilibrium 
expression state S, then —S(t) converges to — S. 

Proof. See Theorem A3 in the Appendix. □ 

Example 3. The gene network given in (17) is used to 
show the validity of Theorem 3. The following two initial 
saturated states 

Si(0) = (11-1111), 

S2(0) = (-1-11-1 -1 -1), 

with S2(0) = — Si(0), were used for dynamics (3). The 
two solution trajectories are found to satisfy Theorem 3. 



Figure 3 gives the projections of the two trajectories onto 
the two-dimensional (Si, Ssj-subspace. 

Corollary 2. In the hypercube [ —1, 1]", the volume of 
the region where points converge to a saturated equilib- 
rium state S is equal to the volume of the region where 
points converge to — S. 

Proof. Considering that S(t) and —S(t) are symmetric 
and have the same distance to the origin (see Figure 3), 
then in [—1,1]" (which is symmetric to the origin) the 
volume of the region where points converge to S is equal 
to the volume of the region where points converge to — S. 
However, if Si and S2 are two saturated fixed point attrac- 
tors, but S2 ^ —Si, for a given W, generally in [ —1, 1]" 
the volume of the region where points converge to Si may 
not be equal to the volume of the region where points 
converge to S2. □ 

Theorem 4. S = 0 is a fixed point, but may not be a 
fixed point attractor for a W having a saturated fixed point 
attractor S. 

Proof. S = 0 is a fixed point because 



Si(t+T) = 





n 


1 -1- exp 


/■=! 



(18) 
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0.8 
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Figure 3 Projections of two trajectories with initial states Si(0) and S2(0)(= -Si(0)), respectively, in Example 3 onto two-dimensional 
(Si,S5)-subspace. 



LetBf = {S e 9^" I ||S|| < e} be an open ball of radius e 
centered at the origin in d]" where e is chosen sufficiently 
small such that there is only a single fixed point 0 within 
The sufficient condition for 0 to be a unique fixed point 
attractor is (see (110) of Theorem Al in the Appendix) 

n 

^|M/y|<2, {i = l,2,...,n). (19) 

/=i 

A saturated fixed point attractor S of W must satisfy (15). 
Then, we have 

n n 

j=i j=i 

= Si ^'/^/ j > /8 > 2, = 1, 2, . . . , «). 

(20) 

Therefore, 0 may not be a fixed point attractor, but an 
unstable fixed point. For an unstable fixed point, there 
exist divergent or both convergent and divergent neigh- 
borhoods of 0. In the convergent region, trajectories will 
be attracted to 0; In the divergent region, trajectories will 
leave from the neighbourhood of 0. Therefore, it is possi- 
ble upon starting from a saturated initial expression state 
S(0) that the solution trajectory arising from the sigmoidal 
function (3) may converge to 0, or first approach 0, but 
then enter the divergent region, and move away from that 
neighborhood. □ 

Theorem 5. A saturated initial expression state S(0) 
converges to a saturated equilibrium expression state S if 



the following condition is satisfied after a finite number k 
of iteration steps of (3) starting from S(0) 

n 

WijSjit >k)> -In [(ffi-D-VCai-D^-l] , i € /+(S), 

;=i 

(21) 

n 

J2 > < In [iai - 1) - - 1)2 - i] , iej- (S), 

(22) 

where 

n 

oii = J2 l^yl (23) 
/=1 

under the constraint that at > 2. 
Proof. See Theorem A4 in the Appendix . □ 

Theorem 5 provides a way to determine all possible 
saturated initial expression states converging to a satu- 
rated equilibrium expression state of a W. The constraint 
Ui > 2 implies that the gene network must contain more 
than one gene if w,y only takes values in [—1, 0, 1]. This is 
always true. 

Example 4. For the gene network given in (17) all 
Ui = 5, and the sufficient condition for a saturated initial 
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expression state converging to a given saturated equilib- 
rium expression state S is then 

XI >k)> 2.0634, / e 7+(S), (24) 

7=1 



/=1 



-2.0634, «e/"(S), 



(25) 



where S(t > k) is the solution of (3) after k iterations 
starting from a saturated state S(0). All 2^ = 64 saturated 
states are used as initial expression states for dynamics (3). 
The results are given in Table 1. 

Using the condition in (24,25), we found that each 
saturated equilibrium expression state for the gene net- 
work given in (17) has 22 saturated initial expression 
states (including the saturated equilibrium expression 
state itself). The remaining 20 saturated initial states, 
which do not satisfy (24,25), converge to 0. 

For saturated initial expression states converging to one 
of the two saturated equilibrium expression states, the 
condition in (24,25) is satisfied starting from the iteration 
number k as either 0 or 1. For saturated initial expression 
states converging to 0, the condition given in (24,25) is 
never satisfied. Therefore, in this gene network, the sig- 
moidal function (3) either does not need to be solved, or 
only needs to be solved just once, to determine which and 
how many saturated initial expression states converge to 
a particular saturated equilibrium expression state. The 
computational effort will be reduced. For a network with 
11 genes, using Theorem 5 has approximately 60% CPU 
time saving compared to completely solving (3). 

Figure 4 gives the projection of a trajectory converging 
to the final state 0 onto the two-dimensional {Si,Ss)- 
subspace. Since 0 may be an unstable fixed point, and the 
values of the elements of S(t) are not exactly zero, but are 
within a small region around 0. Therefore, the trajectory 
is sensitive to computational precision, that determines 
which region the trajectory enters. Continued iteration 
may lead to the trajectory entering the diverging region of 
0 and leaving away from 0. When there are no limit cycles 
for the W, the trajectory must converge to one of the 
two saturated equilibrium expression states. The outcome 
depends on the precision used in the computation. 



Table 1 Number of saturated initial states converging to 
different final states for the gene network given in (1 7) 



Final state 


Number of saturated initial states 


(111-1 -1 -1) 


22 


(-1 -1-111 1) 


22 


(000000) 


20 



Figure 5 shows two trajectories that first approach to 
0 and then leave the neighborhood of 0 and converge 
to one of the two saturated equilibrium expression states 
(111-1 -1 -1) and (-1 -1-111 1), respectively. Such 
a property of 0 for dynamics (3) may not be meaningful 
biologically and causes confusion. This problem will be 
discussed in Section 5. 

3.2 Definition for robustness to noise 

Robustness to noise may be defined as follows for a net- 
work W, for each of its saturated equilibrium expression 
states, and for a specified pair of saturated equilibrium and 
initial expression states, respectively. 

Definition 2. The robustness to noise R„^ of a given 

gene network W is specified as inversely proportional to 
the total number m of fixed points. 



0, if w = 0, 
1/m, otherwise. 



(26) 



The subscript «t denotes 'noise' and 'total'. 



According to the definition, larger values of m corre- 
spond to worse robustness to noise. This definition is 
reasonable because more fixed points a W has, then less 
saturated initial states converging to each of its equilib- 
rium states. Changing an initial state has more chances to 
cause change of the equilibrium state. 

In this definition, m should also include limit cycles. 
Since there is no simple way (like the Banach fixed point 
theorem for the existence of fixed point attractor) to deter- 
mine the existence of limit cycle for sigmoidal functions, 
in the current work, we are unable to include it. The 
robustness i?„^ takes on values in [0,1]. For W given 
in Example 2, there are three fixed points: 0, (1 1 1 — 
1 -1 -1) and (-1 -1-111 1). Therefore, R„, = 1/3, 
which implies that there is a 2/3 probability for having a 
saturated equilibrium expression state change caused by a 
saturated initial expression state change. 

Definition 3. The robustness to noise i?„; of a given sat- 
urated equilibrium expression state S/ for a gene network 
W is specified by the ratio of the number Ni of satu- 
rated initial expression states converging to S;, to the total 
number 2" of possible saturated initial states 



Rn 



■■Ni/2". 



(27) 



The subscripts n and i denote 'noise' and the ith saturated 
equilibrium expression state. For saturated equilibrium 
expression states (111—1 —1 —1) and (—1 —1—111 1) 
in Example 2, R„, = 22/64 « 1/3 (i = 1, 2). 

Ciliberti et al. [13] argued that the pair of saturated equi- 
librium and initial expression states, S(oo) and S(0), play 
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Initial state (-1 -11-11-1) 




40 60 80 

Iteration number, t 



120 



Figure 4 Projection of trajectory onto two-dimensional (Si,S5)-subspace starting from initial state (-1 -11-11 -1) converging to final 
state 0 (Example 4). 



a central role for a viable network, but the variation of ini- 
tial expression state in realistic cases is often mild. They 
define one measure of robustness to noise as the prob- 
ability Rvj that a change in one gene's expression state 
in the saturated initial expression state S(0) leaves the 
unchanged network's saturated equilibrium state S(oo). 
Following this pattern, we also define the measure of 
robustness to noise for a given pair of saturated equilib- 
rium and initial expression states with a viable network as 
follows. 



Definition 4. The robustness to noise of a given 
pair of saturated equilibrium and initial expression states 
S,- and Sj(0) for a viable gene network W is specified by 
the ratio of the number ATy of neighboring saturated initial 
expression states differing from Sy (0) by only one element 
and still converging to S,-, with respect to the total num- 
ber « of possible one element differing saturated initial 
states 



= Nij/n. 



(28) 
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Example 5. Determination of R„ij for the network in 
Example 2. 

The network (17) has two equiUbrium states Si, S2 with 
22 saturated initial states Sy (0) converging to each. The 
Rnij for each pair of S,- and Sy (0) was determined from solv- 
ing dynamics (3) for all n one-element differing saturated 
initial states. The measure R^jS for all possible pairs take 
only two values 1, 1/3 (i.e., Nij only takes value 6 or 2). The 
distribution of R„.^ (i.e., how may pairs have the same J?„,y) 
is given in Table 2. 

For S,- , the seven Sy(0)s in the 22 saturated initial states 
with R„.j = 1 are Sj itself and those differing from S,- by 
only one element; the other fifteen S,(0)s in the 22 satu- 
rated initial states with R„fj = 1/3 are those differing from 
S; by two elements. The Sj(0) differing from S/ by more 
than two elements does not converge to the S,. 

4 Robustness to mutations 

For robustness to noise, we need to find all possible 
saturated equilibrium expression states for a given gene 
network W. In contrast, for robustness to mutations, we 
have the opposite task: for a given saturated equilib- 
rium expression state we need to determine which and 
how many Ws share this saturated equilibrium expression 
state. 

4.1 Conditions under which gene networks share the 
same saturated equilibrium expression state 
Theorem 6. For a given saturated expression state S, all 
possible networks specified by particular Ws having it as a 
saturated equilibrium expression state can be completely 
constructed by solving the following inequalities: 

J2 "^'i- H "^''^ /S, if «e/+(S), (29) 
;6/+(S) 7e/-(S) 

;e/+(S) 7e/-(S) 
under the condition iVy e[ —a, a]. 

Proof. Note that the rows of W are independent. When 

5 is the saturated equilibrium expression state of W, the 
elements of each row (for example, the ith row Wijij = 
1, 2, ... , k)) of W must satisfy either (13) or (14), i.e., 
(29) or (30) which is an inequality with n variables. The 



Table 2 Distribution of R„,j for the network W in Example 2 



Final state 










1 




1/3 


(111-1 -1 -1) 


7 




15 


(-1 -1-111 1) 


7 




15 



inequality is solvable and has an infinite number of solu- 
tions. Those are the desired solutions with each Wij located 
within the required range [ —a, a]. For discrete wtj (only 
taking values [ — 1, 0, 1]), the number of solutions is finite. 
All the solutions can be completely counted and deter- 
mined by solving (29) or (30) . □ 

From (29) or (30), we may draw the following conclu- 
sions: 

1. For / e /"""(S), increasing the value of the first term of 
(29) (if the increase does not make the value of 
larger than the upper bound a) will not violate the 
inequality, i.e., increasing the value of e /"'"(S)) 
will keep the same saturated equilibrium state S. This 
behavior implies that either increasing the activation 
or decreasing the repression influence of active gene / 
on active gene i at equilibrium state S will not change 
the saturated equilibrium state. 

2. For / e /+ (S), decreasing the value of the second term 
of (29) (if the decrease does not make the value of 
smaller than the lower bound —a) will not violate the 
inequality, i.e., decreasing the value of Wij(j e 7~(S)) 
will keep the same saturated equilibrium state S. This 
behavior implies that either decreasing the activation 
or increasing the repression influence of inactive 
gene ) on active gene i at equilibrium state S will not 
change the saturated equilibrium state. 

3. Similarly, for / e /~ (S), either increasing the 
activation or decreasing the repression influence of 
inactive gene j on inactive gene i at equilibrium state 
S will not change the saturated equilibrium state. 

4. For / € /~(S), either decreasing the activation or 
increasing the repression influence of active gene ; 
on inactive gene 2 at equilibrium state S will not 
change the saturated equilibrium state. 

Example 6. Consider the determination of all possible 
gene networks W with the given saturated equilibrium 

expression state 

S = (1 1 1 -1 -1 -1). 

Here, the discrete values [ — 1 , 0, 1 ] are required for the ele- 
ments Wij. Thus, we seek to determine all gene networks 
W sharing the same two saturated equilibrium expression 
states given in Example 2 (due to Theorem 2, — S is also a 
saturated equilibrium expression state). 

First, consider / e /+(S). Set ^ = 5 for the 0.99 
confidence level. In this case, (29) becomes 

3 6 
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Rearrange (31) as 

3 6 
^M',;>5 + ^M',;. (32) 

;=1 /=4 

Note that the second term on the right-hand side of (32) 
takes on integer values from —3 to 3 corresponding to all 
three having the value either —1 or 1, respectively. We 
treat each circumstance separately: 

1. j:wij = -3 

1=4 

In this case, (32) is 

3 

J2wij>2. (33) 

It is easy to see that there is only one choice for 

(W,-4 Wis Wi(,) = (-1 -1 -1) 
and four choices for 

(wn Wi2 Wi3) = (1 1 1), (0 1 1), (1 0 1), (1 1 0). 
Thus, we have four choices for Wy (/ = 1,2, ... ,6) 
wj^ = (1 1 1 -1 -1 -1), 
= (0 1 1 -1 -1 -1), 
= (1 0 1 -1 -1 -1), 
w+ = (1 1 0 -1 -1 -1). 
Here, denotes the Jcth permitted pattern for the 
ith row of W with i e /+(S). 

2. Ewy = -2 

;=4 

In this case, (32) is 

3 

^wy>3. (34) 

It is easy to see that there is only one choice for 

(wa Wi2 w/s) = (111) 

and three choices for 

(W;4 Wis Wi6) = (0 -1 -1), (-1 0 -1), (-1 -1 0). 

Thus, we have another three choices for 
Wij(J = 1, 2, . . . , 6) 

M/+ = (1 1 1 0-1 -1), 

w+ = (1 1 1 -1 0 -1), 
ivt = (1 1 1 -1 -1 0). 



3. T,Wij>-l 
In this case, (32) is 

3 
7=1 



This criterion is impossible because ^ij cannot 
be larger than 3. Therefore, altogether, there are only 
seven choices or permitted rows for Wij when / e /+ (S). 



.,+ 



Ill 


-1 


-1 


-1 


oil 


-1 


-1 


-1 


1 0 1 


-1 


-1 


-1 


110 


-1 


-1 


-1 


111 


0 


-1 


-1 


1 1 1 


-1 


0 


-1 


111 


-1 


-1 


0 



(36) 



Now consider i e /"(S). In this case, (30) becomes 

3 6 

J2 - J2 - = 4. 5, 6). (37) 

Rearrange (37) as 

6 3 
^Wy>5 + ^Wy. (38) 

Note that (38) is the same as (32) except that (w/iw^^ 
vf/s) and iwnWiswi^) interchange their positions. 
Therefore, there are seven choices or permitted rows 
for Wij when i e /~ (S): 



w- 







"-1 


-1 


-1111 






-1 


-1 


-10 11 


M/3 




-1 


-1 


-110 1 






-1 


-1 


-1110 


M/- 




0 


-1 


-1111 






-1 


0 


-1111 


Wj 




-1 


-1 


0 111 



(39) 



(35) 



In the construction of w^, for a given (W;4, w,5, Wjg), all 
possible patterns of (w,i> w/3) are considered. For a 
given {wii,Wi2,Wi3), a similar treatment for {wii,WiStWi() 
was performed. Thus, for any , we can always find a 
differing from it only by a single Wij. This is also true 
for M/^. In this example, and w^(A: = 2, 3, . . . , 7) differ 
from w'^ and only by a single Wij, respectively. 

Since each row of W has seven choices, altogether, 
there are 7^ = 117,649 gene networks, each specified 
by a particular W, sharing the same saturated equilibrium 
expression states 

S = (1 1 1-1-1 -1), 
-S = (-1 -1-1 1 1 1). 

These gene networks sharing the same saturated equi- 
Ubrium expression states are referred to as 'viable' net- 
works (see [20]). This definition is different from that 
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given by Ciliberti et al. [13], where viable networks 
were those sharing a prespecified pair of saturated initial 
and equilibrium expression states. The viable networks 
defined by sharing a pair of saturated initial and equilib- 
rium expression states are a subset of the viable networks 
defined by sharing a saturated equilibrium expression 
state only. 

The analysis here about viable networks implies that 
for a given saturated state, all viable networks having 
it as an equilibrium state must follow certain pat- 
terns, i.e., its rows must be chosen from finite per- 
mitted rows. The permitted rows for a given saturated 
equilibrium state have specific biological meaning and 
reflect the required connectivity patterns of each gene 
to other genes. This restriction distinguishes viable net- 
works for a given equilibrium state from other viable net- 
works with distinct equilibrium states as well as inviable 
networks. 



Suppose that only one can change. From W'^, we see 
that each element Wy in w]*" changing to 0 yields one of 
w^(k = 2,...,7), which is still viable. Other changes are 
not viable. Therefore, the total number of viable single Wy 
changes in w'^ is 

A/'V = 6. (40) 

However, for each in (k — 2, ... ,7) only 0 changing 

to 1 (if; e /+(S)) or -1 (if; € /"(S)) gives w+ which is still 
viable, and other changes are not viable. Thus, the total 
number of viable single w^. changes inw^(k = 2, . . . ,7) is 

N\ = l, ik = 2,...,7). (41) 



4.2 Definitions of robustness to mutation 

The number of viable networks in the example above 
7^ = 117, 649 itself is large, but the total number of pos- 
sible gene networks with « = 6 is 3^ 1.5 x 10^^. The 
fraction of viable gene networks in the total number of 
possible gene networks is 7^/3^^ ^ 7.8383 x 10~^^, even 
smaller than that obtained previously in numerical simu- 
lations [13] for n = 20. Based on the above analysis for 
viable gene networks, it seems plausible to define robust- 
ness to mutation for a given saturated equilibrium state 
as the ratio of the number of viable networks to the total 
number of possible networks for a given n. If so, it would 
appear that none of the viable gene networks is robust to 
mutation because the ratio is very small. 

The latter inference is misleading because most of the 
possible gene networks have no similarity in topology with 
the viable gene networks having a specific saturated equi- 
librium state, and there is rarely a chance that a viable gene 
network will suddenly change to one of them. In normal 
circumstances, the structure perturbations due to muta- 
tion are small and the topology can only change gradually. 
A viable network may experience topology changes step- 
by-step, and in each step, only one W/; changes. Ciliberti 
et al. [13] defined mutational robustness for a viable gene 
network as the fraction of its one-mutant neighbors that 
are also viable, and we follow the same criterion. In the fol- 
lowing discussion, w,; is restricted only to take the discrete 
values [-1,0,1]. 

We will use Example 5 as an illustration. A gene network 
W is viable if and only if its ith row belongs to one of the 
seven rows in W+ (if i e /+(S)) or W' (if i e /"(S)), 
respectively. Since w,; only takes on three values [ — 1, 0, 1], 
each Wij may have two possible changes from its original 
value, and there is a total 2n^ = 2 x 6^ = 72 single Wy 
changes (i.e., 2«'^ = 72 one-mutant neighbors) for any W. 



It can be proved that this is also true for and (k = 
2,..., 7), i.e., 

AT"- =6, N''- = l, {k = 2,...,7). (42) 

We then define robustness to mutation as follows: 

Definition 5. Robustness to mutation for a viable gene 
network W with a specified saturated equilibrium state S; 

is 

0, if W is inviable, 

= \ n (43) 

^ = i:^, if is viable. 

Here the subscripts m and / in denote 'mutation' and 
the j'th saturated equilibrium expression state; and 
N^,. respectively are the total numbers of viable single w,; 
changes (which is also the number of one-mutant viable 
neighbors) of W and its «th row; ATvr is the total number 
of possible single w,; changes of W. For a given W with 
a specified saturated equilibrium state, its robustness to 
mutation can be readily calculated from N"^. of each row. 
In Example 5, since each row of a viable network must be 
one of the seven rows in W+ or W~ , and + = =6, 

N\ = N". = l(k = 2,..., 7), the value of i?^; depends 

on how many and are contained in W. 

Example 7. Some inviable and viable networks Wj^k = 
1, . . . , 8) with respect to the saturated equilibrium expres- 
sion state (111-1 -1 -1) (and (-1 -1-111 1) by 
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Theorem 2) are given below and their and 7?^; ^^re 
given in Table 3. 





- 0 


0 


1 


-1 -1 


-1 - 




- 0 1 


1 


-1 


-1 


-1 - 




1 


0 


1 


— 1 — 1 


— 1 




1 0 


1 


— 1 


— 1 


— 1 


Wi = 


1 


1 


0 


— 1 — 1 


— 1 


W2 = 


1 1 


0 


— 1 


— 1 


— 1 


1 

— i 


1 

— i 


1 

— i 


U 1 


1 
i 


1 1 


1 

— i 


A 

u 


1 
i 


1 
i 




— 1 


— 1 


— 1 


1 C 


1 




— 1 — 1 


— 1 


1 


0 


1 




_ — 1 


— 1 


— 1 


1 1 


0 _ 




_ — 1 — 1 


— 1 


1 


1 


0 _ 




- 1 


1 


1 


-1 -1 


-1 - 




- 1 1 


1 


-1 


-1 


-1 - 




1 


0 


1 


~] ~] 


— 1 




1 0 


1 


— 1 


— 1 


— 1 


W3 = 


1 


1 


0 


— 1 — 1 


— 1 


Wi = 


1 1 


0 


— 1 


— 1 


— 1 


^\ 




_1 


0 1 


I 




_1 


0 


\ 


]^ 




-1 


-1 


-1 


1 0 


1 




-1 -1 


-1 


1 


0 


1 




_ -1 


-1 


-1 


1 1 


0_ 




_ -1 -1 


-1 


1 


1 


1 _ 




- 1 


1 


1 


-1 -1 


-1 - 




- 1 1 


1 


-1 


-1 


-1 - 




1 


1 


1 


-1 -1 


-1 




1 1 


1 


-1 


-1 


-1 


Ws = 


1 


1 


0 


-1 -1 


-1 


Wf, = 


1 1 


0 


-1 


-1 


-1 


-1 


-1 


-1 


0 1 


1 


-1 -1 


-1 


0 


1 


1 




-1 


-1 


-1 


1 0 


1 




-1 -1 


-1 


1 


1 


1 




_ -1 


-1 


-1 


1 1 


1 _ 




_ -1 -1 


-1 


1 


1 


1 . 




- 1 


1 


1 


-1 -1 


-1 - 




- 1 1 


1 


-1 


-1 


-1 - 




1 


1 


1 


-1 -1 


-1 




1 1 


1 


-1 


-1 


-1 


W7 = 


1 


1 


1 


-1 -1 


-1 


Ws = 


1 1 


1 


-1 


-1 


-1 


-1 


-1 


-1 


0 1 


1 


-1 -1 


-1 


1 


1 


1 




-1 


-1 


-1 


1 1 


1 




-1 -1 


-1 


1 


1 


1 




. -1 


-1 


-1 


1 1 


1 _ 




_ -1 -1 


-1 


1 


1 


1 _ 



Wi is inviable because its first row does not belong to 
any row of W+ or W~ . The other s are viable. As men- 
tioned above, the value of Rmt of a viable W depends on 
how many and are contained in the network. From 
W2 to Wg, more and more w'^ and w]" are included. Thus, 
the corresponding Nyrr and Rmf become larger and larger, 
and Ws is the most stable one with N^rr = 36 (i.e., having 
36 one-mutant viable neighbors) and i?^, = 1/2. 

For K = 6 and the saturated equilibrium expression 
states (111-1 -1 -1) and (-1 -1-111 1), the viable 
gene networks can only take seven distinct values of Rm- 
given in Table 3 corresponding to W containing 0,1,. ..,6 
of and wj". Suppose k rows of a viable network are w'^ 
and wj", then the value of its mutational robustness Rmt (k) 
is given by 

6 

Rmtik) = J2 C/(2 X 6^) =[6k + 1(6 - k)] 112 



= {5k + 6) 172, (A: = 0,1,..., 6). 



(44) 



It can be readily proved that the number N{k) of all possi- 
ble viable networks with robustness to mutation equal to 
Rmi{k)is 

N(k) = C^x6^-K (A: = 0,1,..., 6), (45) 
Table 3 The robustness to mutation of Wk 









Wk 




Rm, 




0 


0 


Ws 


21 


21/72 


W2 


6 


6/72 


We 


26 


26/72 


W3 


11 


11/72 


W7 


31 


31/72 


W4 


16 


16/72 




36 


36/72 



where the first term Cg denotes the number of combina- 
tions of k elements taken from six elements, representing 
how many possible positions that k rows of or w]~ can 
take in six rows of W. Each of the remaining 6 — k rows of 
W has six choices from and (k = 2, . . . ,7), and the 
second term 6^~^ gives the total possible number of com- 
binations for the 6 — k rows. Figure 6 gives the distribution 
of Rmr 

Following [13], we also define robustness to mutation, 
Rmij, as follows: 

Definition 6. Robustness to mutation for a viable gene 
network W with specified saturated equilibrium and ini- 
tial states S, and Sj(0) is 



Nw 



AT" 



(46) 



Here, the subscript m and i,j in R„^. respectively denote 
'mutation' and the rth saturated equilibrium expression 
state S; and yth saturated initial expression state (0); N^.. 
is the total numbers of viable single Wy changes of W with 
respect to S/ and Sy(0), which can be obtained by testing 
how many viable networks in share Sy(0). 

Example 8. Determination of Rmij for W2 and W3 in 
Example 7. 

Rmij values were determined for the networks W2 and 

W3 given in Example 7. Both W2 and W3 have the same 
saturated equilibrium states Si, S2. For W2 each S, has 22 
saturated initial states converging to it. For W3 each S; has 
32 saturated initial states converging to it. The numbers 

of one-mutant viable networks sharing the same sat- 
urated equilibrium state S, (z = 1,2) for W2 and W3 are 6 
and 11, respectively (see Table 3). For W2, all 6 one-mutant 
neighbours sharing S, also share all 22 Sy (0)(/ = 1 — 22), 
i.e., Nl,.. =Nlr = 6, R^ij = 6/72 for all 22 pairs. For W^, 

is different not only for distinct S/ but also for distinct 
initial states. For Si, N^..(< ATJ^) takes values: 5,7,11; for 
S2, N^.j (< N^) takes values: 5,6,9,11, respectively. Table 4 
gives the distribution of Rm^. values for W3, i.e., how many 
pairs of S, and Sy(0) take these values. 

4.3 Topology evolution of gene networks 

From the procedure to construct viable networks given in 
Example 6, we know that for each permitted row there 
always exists another permitted row differing by only a 
single Wij from it. Therefore, for a viable network Wi, we 
can always find one or more viable networks, W/s, dif- 
fering by only one from it. The above eight networks 
W^ik = 1, 2, . . . , 8) are an example. They only differ from 
one another as adjacent neighbors with a single changed 
Wij. These changes in topology correspond to the loss of a 
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regulatory interaction {wij — >■ 0), or to the appearance of 
a new regulatory interaction that was previously absent. 
The changes can be represented as a reversible path 

Wi <^ W2 O W3 W4 

^ (47) 
Ws 4^ W7 4^ We 4^ Ws 

In going from W2 to Wi, the gene network no longer 
attains the saturated equilibrium expression state. Thus, 
we may consider Wi as 'dead! In going from W2 to W3, 
however, not only is the saturated equilibrium expres- 
sion state retained, but also the robustness to muta- 
tion becomes higher. Suppose that all possible single 
changes have the same probability, then the gene net- 
work with higher has a greater chance to 'survive'. 
This implies that highly robust topologies can evolve 
from topologies with low robustness through gradual 
Darwinian topological changes or 'natural selection'. 

Ciliberti et al. [13] suggested that all viable networks 
attaining a given gene expression state can be organized 
into a graph whose nodes are networks that differ in their 



Table 4 The distribution oiRmij for the network 



Final state 














5/72 


6/72 


7/72 


9/72 


11/72 


Si 


7 


0 


3 


0 


22 


S2 


1 


6 


0 


3 


22 



topology. Two networks (nodes) in the graph are con- 
nected by an edge if they differ in the value of only one 
regulatory interaction {wij). As proved above, for a viable 
network Wi, we can always find one or more viable net- 
works, WjS differing by only one Wij from it. The number 
of viable neighbors differing by a single Wij for a viable 
network W is simply the value of its ATJ^ (see Table 2). 
Therefore, any two viable networks Wi and Wj with k dif- 
ferent elements Wij can be connected by a path with k 
edges and /c — 1 viable networks between them. For exam- 
ple, W2 and have six different diagonal w^s, and they 
are connected by a path with six edges and five viable 
networks between them. This circumstance implies that 
all viable networks can be organized to comprise a large 
graph which can be easily traversed by a sequence of sin- 
gle Wij changes of network topology. Thus, robustness is 
an evolvable property. To draw this conclusion, a previ- 
ous study performed millions of simulations [13], but the 
analytical treatment here directly leads to this result. 

5 Modified sigmoidal function with threshold 
parameters 

All the results obtained here are based on the sigmoidal 
function model (3) for gene networks. This model is a 
simplified picture, and caution is called for so as to not 
over-interpret the conclusions obtained from our analyt- 
ical treatment. For example, we proved that — S is also 
a saturated equilibrium expression state if S is one; this 
conclusion may not be biologically meaningful. Another 
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conclusion from our analysis is that 0 may be an unsta- 
ble fixed point. Following a previous definition, 0 corre- 
sponds to all genes being 'half-on'. This definition may not 
be appropriate under some circumstances, and instabil- 
ity of 0 introduces difficulty for biological interpretation. 
However, these considerations provide criteria to modify 
the mathematical model, for example, by using the more 
general sigmoidal function proposed in [23] to describe 
network dynamics. To remove — S and 0 from being an 
equilibrium state or a fixed point, the complex sigmoidal 
function given in [23] is unnecessary, we only need to 
slightly modify the sigmoidal model (3) by introducing a 
threshold parameter 0/ [21,22]: 



Si(t + r) = 







1 + exp 








a = 1,2, 


...,«). 



1, 



(48) 



Using (48), 0 is no longer a fixed point because 

2 

Si(t + r) = 







1 -|- exp 


- (z^wijo - e^j 






2 




1 + 


- 1 /= 0, Oi^ 0. 



(49) 



The necessary and sufficient condition for S to be an 
equilibrium state for (48) becomes 



^M/,yS, > p + di, if «e/+(S), 



(50) 



^M/ySy<-jS-h0i, if i&nS). (51) 

Multiplying both sides by —1 and interchanging > and 
< and changing 7+(S) to /"(-S) and /"(S) to /+(-S) 
yield 

n 

7=1 
n 

Y,wii{-Si)> fi-eu if 2e/+(-S). (53) 

7=1 

If6i/ > OforaUj, 

fi-0i<fi + 9i, (54) 
and it is possible that 

n 

Wiii-Sj) tf^+ Oi, if i e /+(-S), (55) 

7=1 



i.e., (50) may not be satisfied, and — S may not be a 
saturated equilibrium state. If Oi < 0 for all i, then 

(56) 



and it is possible that 



'^'7<-'5;) ^ -)S + Bi, if i e r (-S). 



(57) 



7=1 



In this case, (51) may not be satisfied, and — S may not 
be a saturated equilibrium state. 

If Qi < 0 for all i e /+(S) and 0i > 0 for all i e 
/~ (S), then both (50) and (51) may not be satisfied for — S, 
and — S may not be a saturated equilibrium state for W. 
In Example 2, when the model (48) is used with 0i = —1 
(i = 1, 2, 3) and 0, = 2{i = 4, 5, 6), then the network 
given in (17) only has a single saturated equilibrium state 
(111—1—1—1), and all saturated initial states converge 
to it. Thus, the problem reduces to choosing the param- 
eter 9i and giving it biological interpretation. Then, we 
have 

Theorem 7. The necessary and sufficient condition for 

a saturated state S to be an equilibrium expression state or 
a fixed point attractor of the dynamics (48) with a given 
matrix W is 

J2 Wij- J2 "^ij-^i^ p,i{ieJ+(S), (58) 

7e/+(S) 7e/-(S) 

J2 "^v- Yl Wii-Oi<-fi, iiienS), (59) 

/e/+(S) ;e/-(S) 
or equivalently 




Proof. If a saturated state S is a fixed point of (48), it 
must satisfy (58,59), which implies that 



lim Siit) 

t->-00 



1, if 5, = 1 
-1, if5, = -l 



(«' = 1, 2, . . . , n), 
(61) 



i.e., S is a saturated fixed point with 0.99 (/J = 5) or 0.9999 
{p = 10) confidence level. A saturated state S satisfying 
(58,59, or 60) is a fixed point attractor. The proof is given 
in Theorem A5 in the Appendix. □ 

Similarly, we can have 

Theorem 8. For a given saturated expression state S, all 
possible networks specified by particular Ws having it as 
a saturated equilibrium expression state for dynamics (48) 
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can be completely constructed by solving the following 
system of inequalities 

^ M/y - ^ Wij - 6i > p, if / e /+(S), (62) 

;e/+(S) ;£/-(S) 

^ M/y - ^ Wif -6i< -fi, if i e /-(S), (63) 
/■e/+(S) ;e/-(S) 

under the condition e [ —a, a], and Oi e [ —b, b]. 
Proof. The proof is the same as that for Theorem 6. □ 

6 Application to a yeast cell-cycle network 

A simple yeast cell-cycle network shown in Figure 7b with 
11 nodes was proposed by Li et al. [24]. 
The dynamics of the network was defined by Li et al. as 



Si(t + 1) = 



1, ZaijSjit) > 0, 

0, Y:aijSj(t) < 0, 
j 

Si(t), ZaijSjit) = 0, 
j 



(64) 



where 1 and 0 correspond to active and inactive states of 
the gene, i.e., 0 instead of — 1 is used to represent the inac- 
tive state, and fly is Wy in (3). Using model (64), Li et al. 
found that there exist 7 saturated fixed point attractors 
(considering 0 as —1) and all of the 2^^ = 2,048 possible 
saturated initial expression states converge to one of the 
seven fixed point attractors (see Table 5). 

Note that dynamics (64) is different from dynamics (3). 
For "Y^jaijSjit) = 0, dynamics (3) gives Siit +1) = 0, 
not Si{t). Dynamics (3) with the W constructed directly 
from the connectivities in Figure 7b will not give the same 
result as that given by dynamics (64). All the information 
given by the simplified model for yeast cell-cycle network 
(Figure 7b) will be considered as 'available experimental 
information' for budding yeast, and used as an example to 



illustrate our analytical treatment for network construc- 
tion and its robustness analysis. Hereafter, 0 representing 
the inactive state by Li et al. will be replaced by — 1. 
Only the main results are presented here; see the online 
Supplementary information (Additional file 1) for more 
details. 

6.1 Construction of viable networks 

Define the node order from 1 to 11 as specified in Table 5, 
i.e., Cln3 is node 1, MBF is node 2, etc. We first construct 
all viable networks sharing the most stable saturated equi- 
librium expression state, the first fixed point attractor in 
Table 5 

Si = (-1 -1 -1 -1 1 -1 -1 -1 1 -1 -1) 

(65) 

for dynamics (3). According to Theorem 2, these viable 
networks will also share the other saturated equilibrium 
expression state 



-Si. 



(66) 



As mentioned by Li et al. [24], 'the overall dynamic 
properties of the network are not very sensitive to the 
choice of these parameters' (Wjy), but the connectivity pat- 
terns of the network, i.e., the regulatory influence between 
genes (activation, repression, and absence) is important 
for determining gene network robustness. Therefore, we 
restrict vi/,y to only take the discrete values 1 (activa- 
tion), — 1 (repression), and 0 (absence). 

Theorem 6 gives the criterion to construct all of such 
networks W. When only takes values [ — 1, 0, 1], to sat- 
isfy (29,30) each row of W must have five or more nonzero 
elements due to /J > 5. Otherwise, the network would 
not have any saturated equilibrium states. This problem 
occurs not only for networks with less than five genes, 
but also for larger networks with sparse connectivities 




Figure 7 Cell-cycle network of budding yeast (a) and its simplification with only one checkpoint 'cell size' (b). Adapted fronn [24]. 
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Table 5 Fixed point attractors of the cell-cycle network and the number of saturated states (basin size) converging to them 



Basin Cln3 IMBF SBF Cln1,2 Cdhi SwiS Cdc20 ClbS,6 Sid Clb1,2 Mcml 
size 



1,764 


0 


0 


0 


0 


1 


0 


0 


0 


1 


0 


0 


151 


0 


0 


1 


1 


0 


0 


0 


0 


0 


0 


0 


109 


0 


1 


0 


0 


1 


0 


0 


0 


1 


0 


0 


9 


0 


0 


0 


0 


0 


0 


0 


0 


1 


0 


0 


7 


0 


1 


0 


0 


0 


0 


0 


0 


1 


0 


0 


7 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


1 


0 


0 


0 


0 


1 


0 


0 


0 


0 


0 


0 



between genes. For example, Node 1 (Cln3) in Figure 7b is 
a pure 'parent' node, which does not have any regulation 
coming from all other 'children' nodes, i.e., all wy = 0 for 
7 ^ 1, and for Si the condition (30) does not hold: 



or 



;e/+(Si) ;e/-(Si) 



(67) 



To avoid this problem, the factor fi may be introduced 
such that 



W = fiW, 



(68) 



so to satisfy condition (29,30), Wy can only take val- 
ues [—1,0,1] without any restriction on the number of 
nonzero elements in each row of W. For the sake of nota- 
tional simplicity, in the sequel, we still use W instead of 
W, but write dynamics (3) as 



Si(t + x) = a \^WijSj{t) 
2 







1 + exp 









1. (69) 



Conditions (29,30) then become 

J2 "^ij- J2 1, if «e/+(S), (70) 

;e/+(S) jeJ-(S) 

J2 "^'j- J2 "^'f ^ if er(S). (71) 

jef+{S) ;e/-(S) 

For saturated equilibrium state Si, (70,71) may be rewrit- 
ten as 

11 

- J2 - ^ ~ ^'5 ~ ^'9' if = 5' 9' ('72) 
/■=i,;5^5,9 

11 

- X ^ -1 - - w,-9, if J ^5, 9, (73) 



11 



J2 Wij < -1 + Wis + Wi9, if J = 5, 9, (74) 

11 

Wij > 1 + Wis + if « 7^ 5, 9. (75) 

y=i,y^5,9 

Using the condition (74,75), all permitted row pat- 
terns sharing saturated equilibrium state Si for dynamics 
(69) have been completely counted and determined (see 
Additional file 1: Supplementary information). Each row 
has 72,219 permitted patterns. Thus, the total number of 
viable networks sharing saturated equilibrium state Si is 

72,219" 2.7872 x 10^1 

As shown below, for the yeast cell-cycle network, the first 
row of W is restricted to be 

(10000000000), 

then the total number of viable networks for dynamics 
(69) is 

72,219^'' « 3.8594 X 10^^ 

There are many choices of practically relevant networks. 

Similarly, the dynamics with threshold parameters (48) 
is also modified as 



Siit + x) = G \ \^WijSjit) 







1 + exp 


-Pif:wijSj{t)-ei) 







1. 



(76) 



The necessary and sufficient condition to have S as a 
saturated equilibrium state for dynamics (76) becomes 

^ M/y - ^ Wij - di > 1, if i e /+(S), (77) 

;e/+(S) ;e/-(S) 

;e/+(S) ;e/-(S) 



Li and Rabitz EURASIP Journal on Bioinformatics and Systenns Biology 201 4, 201 4:4 
http://bsb.eurasipjournals.eom/content/2014/1/4 



Page 18 of 27 



and for Si (77,78) become 



and 



11 

J2 ^ii + - + ^« + if = 5. 9, (79) 

/■=l,/5^5,9 
11 

Yl, ^ii + ^i> 1 + w<5 + ws, if «' ^ 5, 9. (80) 

/■=l,/5^5,9 

Condition (79,80) can be used to construct all viable net- 
works W for dynamics (76) with a given set of ^^s sharing 
the saturated equilibrium state Si. Since there is no unam- 
biguous biological interpretation for the values of Ou as 
[ — 1, 0, 1], to represent activation, repression, absence for 
Wij, we will not construct all such viable networks here. 

6.2 Construction of yeast cell-cycle networks 

According to the definition for the green and red arrows 
along with the yellow loop in [24], the network directly 
constructed from the connectivities of Figure 7b is 



"-1 0 0 
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0 10 
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0 
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0 0 0 
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0 0 0 
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-1 
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0 0 0 
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1 


-1 



(81) 

Wq, does not satisfy condition (70,71) for any saturated 
state and does not have a saturated equilibrium state for 
dynamics (69). However, Wq will be used as the basis for 
the connectivities of the network to construct networks 
for dynamics (69) with the saturated equilibrium expres- 
sion state Si. The construction of networks reduces to 
satisfying condition (74, 75) or (79, 80) as much as possible 
consistent with experimental observation. 

Two yeast cell-cycle networks for dynamics (69) 



Wi = 
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, (82) 



W2 
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(83) 



have been obtained (the detailed procedure for their con- 
struction can be found in the Additional file 1: Supple- 
mentary information). Wi and W2 differ only for Wg^io. 

We can also use Wq without any change, but introduce 
the threshold parameters 



0 = (di 02 03 04 0s 06 07 0S 09 ^10 ^11 ) 



(84) 



for dynamics (76) satisfying condition (79, 80). One choice 
with the smallest magnitudes of 0,5 



©^=(2 111010000 0) 

is obtained by using 

11 

X! ^i; + = ~1 + ^<5 + W(9. 

;=i,M5,9 
11 

^'i + 0i = 1 + W/5 + Wi9, 

;=i,/5^5,9 



(85) 



ifz = 5,9, (86) 



ifz^5,9. (87) 



6.3 Saturated equilibrium expression states for 
constructed networks 

The saturated equilibrium expression states for a given 

network W in dynamics (69) can be determined by using 
the modified condition of (15) in Theorem 1 




>1, (z = l,2,...,ll). 



(88) 



For n - 11, there are 2" - 2,048 saturated states. All 
of the 2,048 states were tested by condition (88) for Wi 
and W2, respectively, to determine which of them are sat- 
urated equilibrium states for Wi and W2. The test for 
2,048 states took only 0.01 s by Matlab on a Dell Precision 
Workstation T3400. The saturated equilibrium expression 
states for a given network W with threshold vector 0 in 
dynamics (76) can be determined by using the condition 



>1, (« = 1,2,...,11). (89) 
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The saturated equilibrium expression states for Wi, W2, 
and Wo with © are shown below. 

1. Wi 

Wi has two saturated equilibrium expression states 
for dynamics (69) 

51 = (-1 -1 -1 -1 1 -1 -1 -1 1 -1 -1), 

52 = (1 1 1 1-1 1 1 1-1 1 1) 

with 

S2 = -Si. 

2. W2 

Wi has four saturated equilibrium expression states 
for dynamics (69) 

51 = (-1 -1 -1 -1 1 -1 -1 -1 1 -1 -1), 

52 = (1 1 1 1-1 1 1 1-1 1 1), 

53 = (-1 1-1-1 1 ~1 -1 -1 1 -1 -1), 

54 = (1 -1 1 1-1 1 1 1-1 1 1). 

with 

S2 = —Si, S4 = —S3. 

The Si and S3 are just the 1st and 3rd fixed point 
attractors in Table 5. 

3. Wo with 0 given in (85) 

There is only a single saturated equilibrium state for 

dynamics (76) 

Si = (-1 -1-1-1 1-1-1-1 1-1-1). 
6.4 Robustness to noise 

First, the number of saturated initial expression states 
converging to each equilibrium expression state for 
W\, W2 is determined by either directly solving the 
dynamics (69) or using modified condition of Theorem 5 



Table 6 The number of saturated initial states converging 
to different equilibrium states for different gene networks 



P l^XJ mjSjit >k)j>- In [(a,- - 1) - Vfe - 1)2 - 1] , 
ieJ+(S), (90) 

P > ^) j < In [(a; - 1) - Vfe - 1)^ - l] , 

ie/-(S). (91) 



For Wi, the CPU times are 0.8 and 0.3 s, respectively 
to check all 2,048 saturated states, i.e., using Theorem 5 
the CPU time is approximately 41% of that for the direct 
solving of the sigmoidal function. The results are given 
in Table 6. Note that for Wi, W2, no saturated initial 
state converges to the unstable fixed point 0. Therefore, 
in the calculation of we ignore 0 and only consider 



Final state 



Number of saturated initial states 



Wo with 0 



W2 



Si 
S2 
S3 
S4 



2,048 



1,024 
1,024 



979 
979 
45 
45 



the saturated equilibrium states. The resultant robustness 
to noise measures and R„i are given in Table 7. There 
are significant differences between R„. (i = 1, 2, 3, 4) for 
W2. Obviously, the saturated equilibrium states Si,S2 are 
much more stable than S3, S4. 

The robustness to noise i?„,y for each pair of saturated 
equilibrium and initial expression states was calculated. 
The distribution of R„.j, i.e., how many pairs with the same 
value of is given in Tables 8 and 9. 

The results show that Wo with © is completely stable 
for any viable pair; for Wi, there is one neighbour of Sj(p) 
differing at the first element, which causes changes in the 
saturated equilibrium state S;; for W2, the distribution of 
Rmj is divergent, and Si and S2 are much more stable than 
S3 and S4. 



6.5 Robustness to mutation 

The Rmi values have been calculated for Wi, W2, and Wq 
with the 0 given in (85) as shown in Table 10. Note that 
for Si, Rmi is almost the same for Wi, W2, and Wo with 0. 

Robustness to mutation i?^,; for a viable pair of speci- 
fied saturated equilibrium and initial states S; and Sy(0) 
has also been calculated for Wi, W2, and Wq with 0. 
The resultant distribution, i.e., how many pairs having the 
same i?^,^. for W\, W2, and Wq with the 0 given above 
is shown in Figure 8 and Table 512 in Additional file 1: 
Supplementary information. 



7 Conclusion 

Based on the determination of saturated fixed point 
attractors for the sigmoidal function model in (3) with a 
given gene network, W, one can analytically determine 



Table 7 The Robustness to noise R„f and R„f for different 



gene networks 



Network 


^»( 








Si 


S2 


S3 S4 


Wq with G 


1 


1 








1/2 


1/2 


1/2 




W2 


1/4 


0.478 


0.478 


0.022 0.022 
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Table 8 The distribution ofR„y for the networlu with 
© and Wi 



Table 1 0 The robustness to mutation R^i of Wo with 0, 
WuandWz 



Final state 


R„,j (for Wo with 0) 


R„,i (for Wi) 


Final state 


Wo with 0 




Wi 




W2 


11/11 


10/11 






Rmi 


JVV 




Si 


2,048 


1,024 


Si 


137 0.57 


140 


0.58 


143 


0.59 


S2 




1,024 


S2 




140 


0.58 


143 


0.59 








S3 








131 


0.54 








S4 








131 


0.54 



which and how many saturated equilibrium expression 
states exist. Furthermore, for each saturated equilibrium 
expression state of a W, which and how many saturated 
initial expression states converging to it can also be deter- 
mined. These results make it possible to establish the 
robustness of a given gene network to noise without per- 
forming a large number of simulations. Based on the 
necessary and sufficient condition for gene networks to 
share the same saturated equilibrium expression state, one 
can determine all the viable gene networks for a spec- 
ified saturated equilibrium state. This result also makes 
it possible to establish the robustness to mutation for a 
network with a specified saturated equilibrium expression 
state or a specified pair of saturated equilibrium and initial 
expression states. 

The analytical treatment presented here proved that for 
a given saturated state, all viable gene networks having it 
as an equilibrium state must follow certain patterns, i.e., 
the rows of the corresponding W must be chosen from a 
finite number of permitted rows. The permitted rows for 
a given saturated equilibrium state have specific biologi- 
cal meaning and reflect the required connectivity patterns 
of each gene to other genes. This restriction distinguishes 
the viable networks for a given saturated equilibrium state 
from other viable networks with distinct saturated equi- 
librium states as well as inviable networks. The analysis 
also proved, without performing a very large numbers of 
simulations, that all viable networks can be organized as 
a large graph which can be easily traversed by a sequence 
of single Wy changes of network topology. Thus, robust- 
ness is an evolvable property. Highly robust topologies can 
evolve from topologies with low robustness through grad- 
ual Darwinian topological changes or natural selection. 
The analytical treatment presented in this paper may be 
employed not only for robustness analysis but also for the 



Table 9 The distribution of R„ij for the network W2 

Final state Rn^j 

2/11 3/11 4/11 5/11 6/11 7/11 8/11 9/11 10/11 

51 1 12 29 146 791 

52 1 12 29 146 791 

53 1 6 21 7 3 7 

54 1 6 21 7 3 7 



model construction and analysis of other properties for 
gene networks. 

Appendix 

The appendix proves several theorems in the main text. 

Lemma 1 (Banach Fixed Point Theorem [25]). Let 
(X, d) be a non-empty complete metric space with a con- 
traction mapping g : X ^ X. Then g admits a unique 
fixed point x* in X (i.e. g(x*) = x*). Furthermore, x* 
can be found as follows: start with an arbitrary element 
xo e X and define a sequence {x„] by x„ = g(x„_i), then 
x„ X*. 

A map g-.X^Xis called a contraction mapping on X 
if there exists q e[0,l) such that 



d(8(x),g(Y)) < qd(x,Y) 



(92) 



where d denotes the distance, for all x, y e X. 
A continuous function g satisfies the Lipschitz condition 



llg(x)-g(y)||;,<sup||7(t)||p||x-y||^ 
teX 

where /(x) is the Jacobian of g(x). Its ii,j)th entry is 

'dgiix)- 



dxi 



(93) 



(94) 



and ||/(t) \\p is the i^-norm: 
11/11 1 = max ( I /y(t) I 



(=1 




7Mt)/(t)- 



(/•=1,2,--- ,«), (95) 



0 = 1,2,--- ,«), (96) 



(97) 



From (95 to 97) we see that ||/||i is the largest sum of 
the absolute values of the elements in each column; ||/||oo 
is the largest sum of the absolute values of elements in 
each row; and II/II2 is the square root of the largest eigen- 
value for matrix /^(t)/(t). Now define d of g(x) and g(y) 
as the i^-norm of their difference. If g is a contraction 
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mapping, (92) requires that at least one of the L^-norms of 
its Jacobian satisfies 



<q<l. 



(98) 



This condition is sufficient, but not necessary. It is possi- 
ble that one of (95 to 97) is satisfied, but the other two may 
not. Such examples can be constructed. 

Theorem Al. The sufficient condition for the sig- 
moidal function in (3) to have a unique fixed point attrac- 
tor 0 is 



^|M/y|<2, («•= 1,2,...,«). 



(99) 



Proof. Since 0 is a fixed point for any W, to see whether 
it is a unique fixed point attractor, we need to determine 
under what condition (3) is a contraction mapping. 

The (/,7)th entry of the Jacobian (S(/c)) for (3) with 
T = 1 is 

^ 2e-"^W (100) 

1 _|_ 2e~";(*) + e~2"iW '■' 
_ 2 



where 

n 

Ui(k) = ^WijSj(k). 



(101) 
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The matrix form of the Jacobian /'^^+^^(S(^)) is 



Therefore, we have 



Z^+^Hsc^)) = v(k)w 

where V(k) is a diagonal matrix with 
Viiik) = ^ 



The oo-norm is 
r(*+l) 



(102) 



(103) 



||/^'^+^^(S(/c))||oo = max 



i gUi(k) _|_ 2 + e-«iW 

Se[-l,ll« ^ ^ ;=1 

0' = 1,2, ■■■ 



(104) 



Because e"'*,e""'<'') > 0, the maximum of 2/(e"'* +2 + 
g- is given by the minimum of e"' *'^^ + e~"' *^*^^ which 
can be obtained from 



(III:', k ) 

This is true if and only if 



= 0. 



(105) 



(106) 



which yields Ui(k) = 0, e"'^'^^ = e "''^'^^ = 1. The minimum 
for e"'(^' + e-"'(^' is 2, and then 



(107) 



i e"iW _|_ 2 + e-«iW 2' 

Se[-l,ll" 



Figure 9 gives the comparison of + e and 2/(e^ 
2 + fi-^). 



1 

||/«+i>(S(^))||oo = max - ^ I w,y I, a = 1,2,. • • ,«). 

(108) 

To be a contraction mapping requires 

||/'^+i)(S(/c))||oo < 1 (109) 



I.e., 



^|wy|<2, a = 1,2,--- ,«). 

7=1 



(110) 



Notice that 1/2 is the superior value in (107). In many 
cases, not all Si = 0 (i.e., Ui(k)'s are not zero), then the fac- 
tor in (107) has values smaller than 1/2, which implies that 
the condition for a contraction mapping given in (110) 
may be softened as 



^|wy|<2, a = 1,2,--- ,«). 

/=1 



(111) 



The condition in (111) is sufficient, but not necessary. □ 

Theorem A2. A saturated state S satisfying (13) and 
(14) is a fixed point attractor. 

Proof. Suppose S is a fixed point of (3). Let = {S € 
fii" I ||X = S — S|| < e), where e is chosen sufficiently 
small such that there is only a single fixed point S within 
X = 0 is a fixed point in 5^ with representation X 
because S is a fixed point for S in fi^ . If we can prove that 
X = 0 is a fixed point attractor in Bf:, then so is S. 
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Subtracting 5; on the both sides of (3) yields 

Xiit + T) = Si(t + T) - 1 : 





« 


1 + exp 









1-1 







1 + exp 


-f:wij{Sj{t)-i+i) 







1 + exp 



- E (wijXjit) + WijSj) 

7=1 

2 



-2 



-2 



1+e P exp 



2, ze7+(S) 



where jS' > )3. 
The (2,7')th entry of the Jacobian is 



2e ^ exp 


















^1 + e~^' exp 




I)' 









2 "'y 



2e~^' exp 


/■=! 




1 + 2e-^' exp 


/■=! 


+ e"2^' exp 


-2 ZwijXjit) 



e^' exp 



■2 + e-/*' 



e exp 



- EwijXjit) 

7=1 



^M/y, /e/+(S). 



(112) 



(113) 



eP'+2 + e-P' ''' 

Here, the condition X(t) < e » 0 and exp |^E^i ^ij^jit)^ ^ exp j^— ^ 1 were used. The ii-norm of 



row vector /}^^'^^ is 



7=1 



(114) 



When 



eP + 2 + e- 



7=1 



+ 2 + e-^ f 75, if = 5, 



11,014, ifjS = 10, 



(115) 



which is often the case in practice for W satisfying (13) and (14), we have 



||7f+^'(X(0)||i< 1, ieJ+(S). 



(116) 
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The same result can be obtained for i e / (S). Therefore, we have 

||7('+^^(Xa))||oo = max||/f+^\xa))||i < 1, 
i.e., X = 0 is a fixed point attractor, and thus so is S. 



(117) 

□ 



Theorem A3. If an Sit) converges to a saturated equilibrium expression state S, then —Sit) converges to the saturated 
equilibrium expression state — S. 



Proof. For Sit) and —Sit) , the corresponding equations are 
2 

Siit + r) = = =--1, (2 = 1,2,...,«), 



and 





2 






n 




1 + exp 


-Y^WijSjit) 












2 






n 




1 + exp 


-T.wiii-Sjit)) 




/=i 





Siit+x) = 



respectively. It can be proved that 

~Siit+x) = -Siit + x). 
The proof is given below. 

~Siit+x) -- 



1, 0' = 1, 2, ... , n), 







1 + exp 


-llwiji-Sjit)) 







- 1 





n 


1 H- exp 


T.wijSjit) 







1 — exp 


n 

T.wijSjit) 
./=i 


1 + exp 


n 

Y^WijSjit) 


— 2 exp 


n 

J^WijSjit) 

7=1 


1 + exp 


n 

HWySjit) 


1 + exp 


n 

T,WijSjit) 





^1 + exp 


n 

J^WijSjit) 


— 2 exp 


T.WijSjit) 


) 




1 + exp 


n 


) 


^1 + exp 


n 

;=i 


^ ^1 + exp 


n 

- HwijSjit) 

/=! 


) 



1 - 



1 - 



2 ^exp 


T,WijSjit) 


") 




« 








n 




^1 + exp 


J^Wi/Sjit) 


1 + exp 


- T.WijSjit) 


) 





















1 + exp 


- T.^ijSjit) 







-Siit + x). 



(118) 



(119) 



(120) 



(121) 



Li and Rabitz EURASIP Journal on Bioinformatics and Systems Biology 201 4, 201 4:4 
http://bsb.eurasipjournals.com/content/201 4/1 /4 



Page 25 of 27 



Equation (121) implies that there is an one-to-one rela- 
tion between Si(t) and —Si(t). Since S{oo) = S, we have 

S(oo) = -S. (122) 
Therefore, starting from —S(t) will converge to — S. □ 

Theorem A4. A saturated initial expression state S(0) 

converges to a saturated equilibrium expression state S if 
the following condition is satisfied after a finite number k 
of iteration steps of (3) starting from S(0) 

n 

J2 "^'/^/(^ > ^) > - ln[ («i - 1) - - 1)2 - 1] , 
i€j+(S), (123) 



7=1 

i€r(S), (124) 

where 

n 

oli = J2 l^vl (125) 

under the constraint that a/ > 2. 

Proof. First, consider i e /~^{S). If S is the fixed 
point attractor for the initial state S(0), then X = 0 is 
the fixed point attractor for the initial expression state 
X(0) = S(0) - S of (112). This is equivalent to finding 
the condition under which (112) is a contraction mapping. 
According to the Banach fixed point theorem, the suffi- 
cient condition is that the norm of the Jacobian matrix 
satisfies ||/||oo < 1- To prove Theorem A4, we try to seek 
the largest = {X e iW" | ||X|| < a] where a is the upper 
bound to have ||/||oo < 1, i-e., (X(0) || i < 1, for all 

«e/+(S). 

Set 



J = e ^ exp 



/=1 



From (113) the Li-norm of row vector 
represented as 



it+r) 



(126) 



can be 



||/f+^)(X(i))||i 



2y 



(1+7)' 
2a;j 



J2 1'^yl 

<1, «e/+(S), (127) 



which gives the quadratic equation 

y^ + 2(l-ai)y+l = (y-yi)(y-y2)>0. 



where 

yi = (a/ - 1) + ^{at - 1)2 - 1, (129) 
y2 = («/ - 1) - Vfe - 1)2 - 1. (130) 

It is easy to check that yi,y2 are all nonnegative when a/ > 
2, and j2 < 1- 

Equation (128) is valid if and only if y is chosen within 
the two disjoint ranges 

(-00, J2) (yi,00) if > 2, 

(-00, 1) (1, oo) if ai = 2, 

i.e., either smaller than y2 or larger than yi. This implies 
that 



e ^ exp 



or 



7=1 



■- exp 



7=1 



> yi 
< yi 



(131) 



. J < - ln[ {Oil - 1) + - 1)2 - 1] < 0, 
^ -'7^/(^) I , _ _ 1) _ - 1)2 _ 1] > 0. 

(132) 

Using (15), we know that when S(t) is close to the fixed 
point attractor S, 



Y^WijSj>Q, 2e/+(S). 

7=1 

Therefore, we obtain 



(133) 



^ w^Sjif) > - ln[ (at - 1) - ^(ai - 1)2 - 1] . (134) 
7=1 

If S(t = k) satisfies (134), then the trajectory starting from 
S(t = k) will converge to S. Therefore, S{t > k) will also 
satisfy (134). Thus, 

« 

^ mjSjit >k)>- ln[ (ai - 1) - y(ai - 1)2 - 1] . 

7=1 

(135) 

Equation (135) proves the condition for i e /"""(S) given in 
Theorem A4. 
For i e /~(S), (112) becomes 



2 







1 + e^' exp 


- -twijXjit) 




7=1 



(128) 



je/-(S). 
(136) 
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Set 



y = exp 



(137) 



Then the proof procedure is the same as above except that 
fori e /-(S) 



(138) 



when S{t) is close to the fixed point attractor S. Thus, we 
obtain 



J2 ^ijSjit) < - ln[ (ai - 1) + ^{ai - 1)2 - 1] . (139) 



Note that 



- ln[ (a; - 1) + Vioii - 1)2 - 1] 
1 

= ln 



In 



In 



(Ui - 1) + y/(ai -1)2-1 

(a; - 1) - ^{ai -1)2-1 



[(tti - 1) + 7(0-; - 1)2 - 1] [(a,- - l)-y(Q!i- 1)2-1] 
(Q,,. _ 1) _ 7 (a,. -1)2-1 



(Ui - 1)2 - (a; -1)2 + 1 

= ln[ (a; - 1) - Viai - 1)2 - 1] . 
Then we have 



(140) 



J2 ^ijSj(t) < ln[ (at - 1) - - 1)2 - 1] , (141) 
and similarly we obtain 



J2 ^ < {ai-l)-^(ai - 1)2 - 1] , (142) 

i.e., the condition for i e /~ (S) given in Theorem A4. □ 

Theorem A5. A saturated state S satisfying (58) and 
(59) is a fixed point attractor. 

Proof. Suppose S is a fixed point of (48). Let = {S e 
di" I ||X = S — S|| < e ) where 6 is chosen sufficiently small 
such that there is only a single fixed point S within B^. 
X = 0 is a fixed point in 5^ with representation X because 
S is a fixed point for S in 5^ . If we can prove that X = 0 is 
a fixed point attractor in fi^ , then so is S. 



Subtracting on the both sides of (48) yields 

Xiit + T) = Siit + T) - I 

2 

1-1 







1 + exp 









1 + exp 



Y, wijiSjit) -1 + 1) 

2 



1 + exp 



Y,(wijXj(t) + WijSj) ■ 



1 + e ^' exp 



2, ieJ+(S) 



where 



(143) 



(144) 



The following step of proof is exactly the same as that in 
Theorem A2, and will not repeat here. □ 
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